class: inverse,left, middle background-image: url(data:image/png;base64,#background.png) background-size: cover <img src="data:image/png;base64,#LOGO_DIPLOMADO.png" width="500px"/> ##Módulo 2: Análisis espacial y Percepción Remota satelital óptica ### Análisis espacial MatÃas Olea <br> <a href="https://orcid.org/0000-0003-0194-7784"> ORCID </a><br> matias.olea@pucv.cl</a><br> .large[<b><a href="https://www.pucv.cl/uuaa/site/edic/base/port/labgrs.html">LabGRS</a> | Septiembre 2024</b>] <br> --- class: center,middle background-image: url(data:image/png;base64,#labgrs_logo.png) background-size: 35% --- ## Contenidos .pull-left[ 1) Operaciones de consulta a un ráster. 2) Operaciones de transformación a un ráster. 3) Operaciones de alteración a un raster. 4) Operaciones de extracción de información. 5) Análisis espacial entre ráster. ] .pull-right[ <img src="data:image/png;base64,#https://raw.githubusercontent.com/allisonhorst/stats-illustrations/main/rstats-artwork/r_rollercoaster.png" width="650px"/> ] --- ## Modelo ráster <center> <img src="data:image/png;base64,#modelo_raster1.png" width="900px"/> </center> --- ## Formatos ráster Los modelos ráster pueden ser almacenados en múltiples formatos, los más comunes son: -- - **GeoTIFF** (Geographical Tagged Image File Format). Fuente de datos en formato Tiff, con la particularidad de que es <u>capaz de almacenar información de referencia espacial</u>. -- - **GMLJP2** (Geographic Markup Language JPEG2000). Muy parecido al formato JPEG pero con <u>mayores beneficios a la hora de comprensión de los datos</u> sin perder información o calidad. -- - **NetCDF** (Network Common Data Form). Fuente de datos creada para el almacenamiento de información cientÃfica en varias dimensiones. Esta fuente de información se asocia principalmente a datos de tipos climáticos o productos satelitales. Su <u>gran capacidad de información</u> la hacen muy útil para compartir grandes volúmenes de información. -- - **HDF** (Hierarchical Data Format). Al igual que el NetCDF, son fuentes creadas para almacenar <u>grandes volúmenes de información</u> cientÃfica espacial. -- - **HGT** Formato asociado a modelo de elevación SRTM (Shuttle Radar Topography Mission Height) <u>caracterizado por almacenar valores enteros positivos y negativos de 2 bits</u>. --- ## Operaciones de consulta en R -- ``` r library(terra) elevacion_3672 <- rast("C:/TuDirectorio/") plot(elevacion_3672) ``` <img src="data:image/png;base64,#DIPGEOPR_02_1_files/figure-html/unnamed-chunk-2-1.png" width="100%" /> --- ## Operaciones de consulta en R **Clase de objeto** ``` r class(elevacion_3672) ``` ``` ## [1] "SpatRaster" ## attr(,"package") ## [1] "terra" ``` -- **Dimensiones** ``` r c(nrow(elevacion_3672), ncol(elevacion_3672), nlyr(elevacion_3672)) ``` ``` ## [1] 3601 3601 1 ``` -- **Resolución** ``` r res(elevacion_3672) ``` ``` ## [1] 0.0002777778 0.0002777778 ``` --- ## Operaciones de consulta en R **Extensión** ``` r ext(elevacion_3672) ``` ``` ## xmin xmax ymin ymax ## -72.00014 -70.99986 -36.00014 -34.99986 ``` --- ## Operaciones de consulta en R **Sistema de referencia y Coordenadas** ``` r crs(elevacion_3672) ``` ``` ## [[1]] ## [1] "GEOGCRS[\"WGS 84\"," ## [2] " DATUM[\"World Geodetic System 1984\"," ## [3] " ELLIPSOID[\"WGS 84\",6378137,298.257223563," ## [4] " LENGTHUNIT[\"metre\",1]]]," ## [5] " PRIMEM[\"Greenwich\",0," ## [6] " ANGLEUNIT[\"degree\",0.0174532925199433]]," ## [7] " CS[ellipsoidal,2]," ## [8] " AXIS[\"geodetic latitude (Lat)\",north," ## [9] " ORDER[1]," ## [10] " ANGLEUNIT[\"degree\",0.0174532925199433]]," ## [11] " AXIS[\"geodetic longitude (Lon)\",east," ## [12] " ORDER[2]," ## [13] " ANGLEUNIT[\"degree\",0.0174532925199433]]," ## [14] " ID[\"EPSG\",4326]]" ``` --- ## Operaciones de consulta en R **Tipo de dato** ``` r datatype(elevacion_3672) ``` ``` ## [1] "INT2S" ``` **_Codificación_**
--- <center><img src="data:image/png;base64,#datatype.png" height="600px" /></center> --- <center><img src="data:image/png;base64,#res_datatype.png" height="600px" /></center> --- ## Operaciones de consulta en R Ahora, esta información que consultamos una por una, la podemos obtener de forma fácil y resumida solo imprimiento nuestro objeto: ``` r elevacion_3672 ``` ``` ## class : SpatRaster ## dimensions : 3601, 3601, 1 (nrow, ncol, nlyr) ## resolution : 0.0002777778, 0.0002777778 (x, y) ## extent : -72.00014, -70.99986, -36.00014, -34.99986 (xmin, xmax, ymin, ymax) ## coord. ref. : lon/lat WGS 84 (EPSG:4326) ## source : s36w072.hgt ## name : s36w072 ``` --- ## Operaciones de consulta en R **Consulta booleana** ``` r consulta_250 <- elevacion_3672 values(consulta_250) <- as.numeric(values(consulta_250) <= 250) plot(consulta_250) ``` <img src="data:image/png;base64,#DIPGEOPR_02_1_files/figure-html/unnamed-chunk-16-1.png" width="100%" /> --- ## Operaciones de alteración ** Uniones y Mosaicos** Las uniones y mosaicos nos permiten unir más de un geodato ráster. Por principio, la union/merge nos permite unir cuando no hay zonas de traslape entre los ráster, y mosaico une cuando hay zonas de traslape. Teniendo en cosideración ello, la función **mosaic()** de terra, es más rápida que **merge()**. <center><img src="data:image/png;base64,#mosaico.png" width="600px"/></center> --- ## Operaciones de alteración ** Uniones y Mosaicos** ``` r coleccion_elev <- list.files("TuDirectorio/", glob2rx("*hgt"), full.names=T) %>% sprc() mosaico_elev <- mosaic(coleccion_elev) plot(mosaico_elev) ``` <img src="data:image/png;base64,#DIPGEOPR_02_1_files/figure-html/unnamed-chunk-18-1.png" width="100%" /> --- ## Operaciones de alteración **Recorte, Clip o Crop** El recorte signica sustraer parte de la matriz ráster según una extensión definida. Podemos crear una extensión ingresando las <u>coordenas lÃmite</u> (xmin, xmax, ymin, ymax), usando la <u>extensión de un geodato</u>, o <u>dibujando una extensión</u> sobre un plot usando la función **draw()**. <center><img src="data:image/png;base64,#clip.png" width="600px"/></center> --- ## Operaciones de alteración **Recorte, Clip o Crop** ``` r cuenca_maule <- vect("TuDirectorio/Cuenca_Rio_Maule.gpkg") %>% project(crs(mosaico_elev)) recorte_elev <- mosaico_elev %>% crop(ext(cuenca_maule)) plot(recorte_elev) plot(cuenca_maule, add=T) ``` <img src="data:image/png;base64,#DIPGEOPR_02_1_files/figure-html/unnamed-chunk-20-1.png" width="100%" /> --- ## Operaciones de alteración **Máscara o Mask** La aplicación de una máscara o enmascaramiento implica el sustraer los valores de un ráster dentro de los lÃmites de un geodato, reclasificanto los pixeles fuera de estos lÃmites como "No Dato", "not available" o "NA". Se aconseja recortar y enmascarar al área de interes. <center><img src="data:image/png;base64,#mascara.png" width="900px"/></center> --- ## Operaciones de alteración **Máscara o Mask** Un elemento importante a considerar es hacerse la pregunta, ¿Qué ocurre cuando el lÃmite de mi área de interés toca parte del pixel pero no lo cubre por completo? ¿lo incluye o no? Existe un parámetro de tipo lógico dentro la función **mask()** llamado "touches". Cuando touches es **T** si toca el pixel, queda incluido en la máscara, si touches es **F** el pixel debe estar contenido en un 100% dentro del área de interes para ser incluido. <center><img src="data:image/png;base64,#touches.png" width="400px"/></center> --- ## Operaciones de alteración **Máscara o Mask** ``` r mascara_completa <- mosaico_elev %>% mask(cuenca_maule) mascara_crop <- recorte_elev %>% mask(cuenca_maule) ``` <img src="data:image/png;base64,#DIPGEOPR_02_1_files/figure-html/unnamed-chunk-22-1.png" width="100%" /> --- ## Operaciones de alteración **Reclasificación** La reclasificación significa reasignar nuevos valores a las celdas del ráster. Por lo general se utiliza cuando queremos pasar de valores continuos a valores discretos para agrupar por rangos en categorÃas especÃficas. <center><img src="data:image/png;base64,#reclass.png" width="900px"/></center> --- ## Operaciones de alteración **Reclasificación** ``` r # Opción 1 reclass_1 <- mascara_crop # guardamos en otro objeto por si nos equivocamos reclass_1[reclass_1 <= 500] <- 0 reclass_1[reclass_1 > 3000] <- -4 reclass_1[reclass_1 > 2000] <- -3 reclass_1[reclass_1 > 1000] <- -2 reclass_1[reclass_1 > 500] <- -1 ``` <img src="data:image/png;base64,#DIPGEOPR_02_1_files/figure-html/unnamed-chunk-24-1.png" width="100%" /> --- ## Operaciones de alteración **Reclasificación** ``` r # Opción 2 intervalos <- c(-Inf,500,0,500,1000,1,1000,2000,2,2000,3000,3,3000,Inf,4) matriz_reclass <- matrix(intervalos, nrow=5, ncol=3, byrow=T) ``` ``` ## [,1] [,2] [,3] ## [1,] -Inf 500 0 ## [2,] 500 1000 1 ## [3,] 1000 2000 2 ## [4,] 2000 3000 3 ## [5,] 3000 Inf 4 ``` --- ## Operaciones de alteración **Reclasificación** ``` r # Opción 2 reclass_2 <- classify(mascara_crop, matriz_reclass) ``` <img src="data:image/png;base64,#DIPGEOPR_02_1_files/figure-html/unnamed-chunk-28-1.png" width="100%" /> --- ## Operaciones de alteración **Remuestreo** El remuestreo es una técnica que se utiliza para transformar la resolución de nuestro dato. Cuando queremos reducir el tamaño del pixel original, estamos hablando de **downscaling**, y cuando queremos aumentar el tamaño, **upscaling**. Comunmente recurrimos al remuestreo cuando tenemos fuentes de <u>información con diferente resolución nativa</u>. Por ejemplo la información climática grillada del CR2 (CR2MET) tiene un pixel de 5 km. El modelo de elevación que estamos usando tiene una resolución de 30m. Las imagenes del sensor MODIS tienen una resolución de 250 m en algunos productos. Por otra parte, tambien utilizamos el remuestreo cuando se realizan reproyecciones. Tópico que veremos en las operaciones de transformación. Entre los métodos de remuestro más comunes están: Vecino más cercano (Nearest neighbor), Interpolación bilineal (Bilinear), Convolución cúbica (Cubic convolution) y por mayorÃa (Majority). --- ## Operaciones de alteración **Remuestreo** .pull-left[ **Vecino mas cercano (Nearest Neighbor)** Este método asiga como valor de salida, el valor de la celda original cuyo centro esté posicionado en la nueva celda. Este método se utiliza cuando no queremos alterar el valor original de entrada sino que solo reposicionarlo, es el caso de las variables discretas, como por ejemplo los usos de suelo. ] .pull-right[ <center><img src="data:image/png;base64,#vecino.png" width="500px"/></center> ] --- ## Operaciones de alteración **Remuestreo** .pull-left[ **Interpolacion bilinear (Bilinear)** Este método no solo toma en cuenta el pixel de origen que ocupa el centro del nuevo pixel, sino que toma en consideración otros 3 y promedia los valores en porcentajes equivalentes a la superficie que cubren del nuevo pixel. Este método tiene a suavizar las transiciones de valores, se utiliza a menudo cuando hay ruido en los datos de entrada. ] .pull-right[ <center><img src="data:image/png;base64,#BILINEAR.png" width="500px"/></center> ] --- ## Operaciones de alteración **Remuestreo** .pull-left[ **Convolución Cúbica (Cubic Convolution)** Este método es parecido al bilineal pero en lugar de considerar 4 pixeles, considera 16. Es utilizado cuando existe exceso de ruido y tiene como resultado un suavizado mayor que el bilineal. ] .pull-right[ <center><img src="data:image/png;base64,#CUBICA.png" width="500px"/></center> ] --- ## Operaciones de alteración **Remuestreo** .pull-left[ **Mayoria (Majority)** Este método es muy similar al vecino más cercano, la gran diferencia es que en lugar de adquirir el valor del pixel que se posiciona en el centro del nuevo, adquiere el valor del pixel que cubre mayor parte de la superficie del nuevo. ] .pull-right[ <center><img src="data:image/png;base64,#MAYORIA.png" width="500px"/></center> ] --- ## Operaciones de alteración **Remuestreo** <center><img src="data:image/png;base64,#plot_dem.png" width="1000px"/></center> --- ## Operaciones de alteración **Remuestreo** <center><img src="data:image/png;base64,#perfil.png" width="1000px"/></center> --- ## Operaciones de alteración **Remuestreo** <center><img src="data:image/png;base64,#plot_lc.png" width="1000px"/></center> --- ## Operaciones de transformación **Reproyección** La reproyección implica reposicionar en el espacio XY nuestro geodato teniendo en consideración otro punto de referencia, proyección y/o coordenadas. Si no conocemos el sistema de referencias y coordenadas, podemos revisar en [EPSG.IO](https://epsg.io/) <center><img src="data:image/png;base64,#epsg.png" width="600px"/></center> --- ## Operaciones de transformación **Reproyección** <center><img src="data:image/png;base64,#Resampling-methods.png" width="550px"/></center> .footnote[ <span style="font-size:9pt"> Fuente: T. Outtara et al (2004); Remote Sensing and Geosciences </span> ] --- ## Operaciones de transformación **Reproyección** ``` r reclass_utm <- project(reclass_2, "EPSG:32719", method="near", res=30) crs(reclass_utm, describe=T) ``` ``` ## name authority code ## 1 WGS 84 / UTM zone 19S EPSG 32719 ## area ## 1 Between 72°W and 66°W, southern hemisphere between 80°S and equator, onshore and offshore. Argentina. Bolivia. Brazil. Chile. Colombia. Peru ## extent ## 1 -72, -66, -80, 0 ``` --- ## Operaciones de transformación **Cambio de tipo de dato** ``` r datatype(elevacion_3672) ``` ``` ## [1] "INT2S" ``` ``` r summary(elevacion_3672) ``` ``` ## s36w072 ## Min. : 3.0 ## 1st Qu.: 141.0 ## Median : 240.0 ## Mean : 403.1 ## 3rd Qu.: 468.0 ## Max. :3031.0 ``` ``` r writeRaster(elevacion_3672, filename="TuDirectorio/Elevacion.tif", datatype="INT1U") ``` --- ## Operaciones de transformación **Vectorización: poligonización** ``` r reclass_pol <- as.polygons(reclass_utm, dissolve=T, values=T) ``` <img src="data:image/png;base64,#DIPGEOPR_02_1_files/figure-html/unnamed-chunk-33-1.png" width="100%" /> --- ## Operaciones de transformación **Vectorización: contornos o isolineas** ``` r nivel <- c(250,500,1500,2500,3500) elev_contorno <- as.contour(mascara_crop, levels=nivel) # creamos vector de lineas plot(mascara_crop) terra::contour(mascara_crop, levels=nivel, add=T) # nos permite plotear encima con etiqueta ``` <img src="data:image/png;base64,#DIPGEOPR_02_1_files/figure-html/unnamed-chunk-35-1.png" width="100%" /> --- ## Operaciones de extracción de información. **Cálculo de área** ``` r frec_reclass <- freq(reclass_utm) frec_reclass$area <- (frec_reclass$count * 900)/1000000 ```
--- ## Operaciones de extracción de información. **Cálculo de área** .pull-left[ Abra el Landcover de la cuenca del Maule [(Zhao et al., 2016)](https://doi.org/10.1016/j.rse.2016.05.016), y calcule el área para cada una de las coberturas de uso de suelo descritas. Utilice como referencia la tabla adjunta que detalla el significado de los valores del producto. - ¿Cuál es la clase con mayor superficie? - ¿Cuál es la superficie más pequeña en km2? ] .pull-right[ <center><img src="data:image/png;base64,#landcover.png" width="350px"/></center> ] --- ## Operaciones de extracción de información. **Cálculo de área** ``` r landcover_maule <- rast("TuDirectorio/Landcover_Cuenca_Maule.tif") frec_landcover <- freq(landcover_maule) frec_landcover$value[which.max(frec_landcover$count)] # valor 410 min(frec_landcover$count)*900/1000000 # 0.0225 km2 ``` ``` ## [1] 410 ``` ``` ## [1] 0.0225 ``` --- ## Operaciones de extracción de información. **Extracción de valores de los pixeles** ``` r hist(mascara_crop, main="Histograma elevación cuenca Maule") ``` <img src="data:image/png;base64,#DIPGEOPR_02_1_files/figure-html/unnamed-chunk-40-1.png" width="100%" /> --- ## Operaciones de extracción de información. **Extracción de valores de los pixeles** ``` r valores_elev <- values(mascara_crop) %>% na.omit # omitimos las celdas NA mean(valores_elev) ``` ``` ## [1] 875.3846 ``` ``` r median(valores_elev) ``` ``` ## [1] 388 ``` --- ## Operaciones de extracción de información. **EstadÃstica zonal** Para este ejercicio extraeremos el valor promedio de elevación por ciudad de la provincia de Talca de la región del Maule. ``` r comunas_maule <- vect("TuDirectorio/Region_Maule_Difrol_latlong.gpkg") provincia_talca <- subset(ciudades_maule, ciudades_maule$NOM_PROVIN == "TALCA") media_elev_comunas <- terra::extract(mosaico_elev, provincia_talca, fun=mean, method="simple", touches=T) media_elev_comunas$ciudad <- provincia_talca$URBANO ``` --- ## Operaciones de extracción de información. **EstadÃstica zonal**
--- ## Análisis espacial entre ráster. **Algebra de mapas o Algebra ráster** <center><img src="data:image/png;base64,#algebra_mapas.png" width="800px"/></center> --- ## Análisis espacial entre ráster. **Algebra de mapas o Algebra ráster** En un caso hipotetico donde haya que definir sitios de interés cientifico o para la conservación, se definieron las variables de EVI (Indice de Vegetación mejorada), la elevación y el Landcover de la cuenca del rÃo Maule para generar zonas de interés. Teniendo en consideración un puntaje de 0 para zonas menos aptas, y 4 para más aptas, reclasifique sus variables de acuerdo a los siguientes criterios: <center><img src="data:image/png;base64,#matriz_emc.png" width="600px"/></center> Una vez reclasificado, multiplique las 3 variables y revise cuales son las áreas con mayor puntaje. --- ## Análisis espacial entre ráster. **Algebra de mapas o Algebra ráster** ``` r ## Apertura de Datos elev_maule <- rast("TuDirectorio/ELEV_CUENCA_MAULE.tif") land_maule <- rast("TuDirectorio/Landcover_Cuenca_Maule.tif") evi_maule <- rast("TuDirectorio/EVI_CUENCA_MAULE.tif") ## Matrices de reclasificación class_elev <- c(-Inf,0,0,0,1000,1,1000,2000,2,2000,3000,3,3000,Inf,4) matriz_elev <- matrix(class_elev, nrow=5, ncol=3, byrow=T) class_evi <- c(-Inf,0,0,0,0.15,1,0.15,0.3,2,0.3,0.45,3,0.45,Inf,4) matriz_evi <- matrix(class_evi, nrow=5, ncol=3, byrow=T) class_land <- c(0,150,0,200,239,4,240,260,1,300,390,2,400,480,4,500,700,3,799,Inf,0) matriz_land <- matrix(class_land, nrow=7, ncol=3, byrow=T) ``` --- ## Análisis espacial entre ráster. **Algebra de mapas o Algebra ráster** <img src="data:image/png;base64,#DIPGEOPR_02_1_files/figure-html/unnamed-chunk-48-1.png" width="100%" /> --- ## Análisis espacial entre ráster. **Algebra de mapas o Algebra ráster** .pull-left[ Boxplot de los pixeles ``` r boxplot(Area_interes,horizontal=F) ``` <img src="data:image/png;base64,#DIPGEOPR_02_1_files/figure-html/unnamed-chunk-49-1.png" width="100%" /> ] .pull-right[ Cálculo del percentil 95 ``` r quantile(values(Area_interes) %>% na.omit(), probs = 0.95) ``` ``` ## 95% ## 24 ``` ] --- ## Análisis espacial entre ráster. **Algebra de mapas o Algebra ráster** <img src="data:image/png;base64,#DIPGEOPR_02_1_files/figure-html/unnamed-chunk-51-1.png" width="100%" /> --- ### BibliografÃa 2022 Jerry Davis. Introduction to Environmental Data Science. Capitulo 8.Raster Spatial Analysis. https://bookdown.org/igisc/EnvDataSci/raster-spatial-analysis.html 2022 Robert J. Hijmans. terra.pdf. https://cran.r-project.org/web/packages/terra/terra.pdf 2012 Essentials of Geographic Information Systems. Chapter 8.Geospatial Analysis II: Raster Data. https://saylordotorg.github.io/text_essentials-of-geographic-information-systems/s12-geospatial-analysis-ii-raster-.html --- class: inverse middle 